classdef RootHelper< handle
    %UNTITLED12 此处显示有关此类的摘要
    %   此处显示详细说明

    properties (Constant)
        ww = 2;
        wh = 2;

        trueW = 240 - Helper.ww;
        trueH = 150 - Helper.wh;
        whRatio = Helper.trueW / Helper.trueH;

        args = {};
    end
    methods(Static = true)
        function [BWWs, data, tt] = ExtractRoots(img8, args)
%             addpath('F:\UAV\马铃薯');
            bShow = args.bShow;
            % % %             rgb = double(img8);
            % % %             red = rgb(:, :, 1);
            % % %             green = rgb(:, :, 2);
            % % %             blue = rgb(:, :, 3);
            % % %             gray = rgb2gray(img8);
            % % %             img = rescale(gray);%%
            
            lab = double(rgb2lab(img8));
            lll = lab(:, :, 1);
            aaa = lab(:, :, 2);
            bbb = lab(:, :, 3);
            img = rescale(lll);%%257-2.jpg 处理好
            aaa = rescale(aaa);
            bbb = rescale(bbb);

            % % % %             ycr = double(rgb2ycbcr(img8));
            % % % %             yyy = ycr(:, :, 1);
            % % % %             ccc = ycr(:, :, 2);
            % % % %             rrr = ycr(:, :, 3);           
            % % % %             ccc = rescale(ccc);
            % % % %             rrr = rescale(rrr);
            %             img = rescale(yyy);%%257-2.jpg 处理不好


            % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
            BW = RootHelper.ExtractWholeMask(img, false);%%Step 1:获得全部种子区域
            % % % % % % % % % % % % % %
% %             BW = RootHelper.SetConnectRegion(BW, bShow);%%Step 2a:  对弱连通区域 加强连接 2022-02-09好像不再需要，需2验证
            BWWs = RootHelper.GetSeedsMask(BW, bShow);%%Step 2b:  提取各个种子根系的单独 掩膜图像


            %%%%%%图像特征分析！ OK
            % % % % % % % % % % % % % % % % % OK, 单独显示 ccc 与 rrr分量
%             bSave = true;
%             bSave = false;     
            
% % % % % %             DetailHelper.saveImage(ccc, bSave, args.fileSaveCCI)%%把aaa，bbb，ccc，rrr的灰度图画出
% % % % % %             DetailHelper.saveImage(rrr, bSave, args.fileSaveRRI)
% % % % % %             DetailHelper.saveImage(aaa, bSave, args.fileSaveAAI)
% % % % % %             DetailHelper.saveImage(bbb, bSave, args.fileSaveBBI)
% % % % %             return
            
% %             bSave = true;
% % % %             bSave = false;             
% %             DetailHelper.showHistgram(ccc, BW, '直方图分割 ccc', bSave, args.fileSaveCCH)%%把直方图画出
% %             DetailHelper.showHistgram(rrr, BW, '直方图分割 rrr', bSave, args.fileSaveRRH)
% %             DetailHelper.showHistgram(aaa, BW, '直方图分割 aaa', bSave, args.fileSaveAAH)
% %             DetailHelper.showHistgram(bbb, BW, '直方图分割 bbb', bSave, args.fileSaveBBH)            
% % %             return 
            
% %             bSave = true;
% %             bSave = false;     
% %             DetailHelper.showTwoThresh(ccc, BW, '二阈值分割ccc', bSave, args.fileSaveCC2)%%把二阈值分割画出
% %             DetailHelper.showTwoThresh(rrr, BW, '二阈值分割rrr', bSave, args.fileSaveRR2)
% %             DetailHelper.showTwoThresh(aaa, BW, '二阈值分割aaa', bSave, args.fileSaveAA2)
% %             DetailHelper.showTwoThresh(bbb, BW, '二阈值分割bbb', bSave, args.fileSaveBB2)            
% %             return
            
            %%%%%%%%%%三阈值
% % % % %             bSave = true;
% % % % % % % %             bSave = false;
% % % % %             DetailHelper.showThreeThresh(ccc, BW, '三阈值分割 ccc', bSave, args.fileSaveCCC)%%把三阈值分割画出
% % % % %             DetailHelper.showThreeThresh(rrr, BW, '三阈值分割 rrr', bSave, args.fileSaveRRR)
% % % % %             DetailHelper.showThreeThresh(aaa, BW, '三阈值分割 aaa', bSave, args.fileSaveAAA)
% % % % %             DetailHelper.showThreeThresh(bbb, BW, '三阈值分割 bbb', bSave, args.fileSaveBBB)            
% % % % %             return
            
            % % % % % % % % % % % % % % % % % % % % % % 聚类 分析
            % % %             DetailHelper.kmeanColors(img8);%%%用RGB进行聚类，没有用
            % % %
            % % %             return;            
            % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %       
            % % %             BWWs = RootHelper.RegineEachSeedRegion(BW, BWWs, ccc, rrr);%%Step2c:没用
            % %             return;
            % % % % % % % % % % % % % % % % % % OK            
%             BWWs = RootHelper.getColorFromBox(BWWs, ccc, rrr, img);%%Step 3a:%%截取不同颜色分量，OK
            BWWs = RootHelper.getColorFromBox(BWWs, aaa, bbb, img);%%Step 3a:%%截取不同颜色分量，2022-05-25 尝试，aaa里面有三阈值结构
%             return;
                        
            threshC3 = multithresh(aaa(BW), 2);%%三阈值 2022-05-24 添加
            threshR3 = multithresh(bbb(BW), 2);%%三阈值 2022-05-24 添加
            threshC2 = multithresh(aaa(BW), 1);%%二阈值 2022-05-25 添加
            threshR2 = multithresh(bbb(BW), 1);%%二阈值 2022-05-25 添加
% %             [threshC3, threshR3, threshC2, threshR2, threshNoR] = DetailHelper.GetMultiThresh(aaa, bbb, BW);
            BWWs.threshC3 = threshC3;%%返回是个数组， threshC(1)<threshC(2)
            BWWs.threshR3 = threshR3;
            BWWs.threshC2 = threshC2;
            BWWs.threshR2 = threshR2;
% % %             BWWs.threshNoR = threshNoR;%% threshNoR没用
            
            BWWs = RootHelper.getSeedRegion(BWWs, bShow);%%Step 3b:进行阈值分割，获得多个种子区域的初始部分2022-04-04
           
%             return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            BWWs = RootHelper.RefineSeedRegion(BWWs, bShow);%%Step 3c:对连通的两个种子区域进行分割，获取不相连的种子区域
%              return;
            BWWs = RootHelper.getBudAndRoot(BWWs, bShow, args);%%Step 3d:对每个种子区域  进行各个部分的分割2022-04-04
% %            return;
    
            BWWs = RootHelper.refineFalseSeed(BWWs, args);%%Step 3e:去除 冗根冗芽，多余的2022-04-04
 % %             return; %%可以从Step 3e中绘制 Seed-root-bud的直方图         
            args.img8 = img8;
            BWWs = RootHelper.reNumberSeeds(BWWs, bShow, args);%%Step 3f:对种子重排顺序 2022-04-04
% % %             return; 
            
            %%再保存一次最后的直方图            
% % %             RootHelper.calcPixelsEachRegion(BWWs, args);    
% % %             return; %%可以从Step 3e中绘制 Seed-root-bud的直方图      

            % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %OK 具体去除每一部分
            BWWs = RootHelper.getBudAndSeed(BWWs, bShow);%%Step 4a
%             return;
            
            BWWs = RootHelper.getRefinedBuds(BWWs, bShow);%%Step 4b
%             return;
            BWWs = RootHelper.getMajorRoots(BWWs, bShow);%%Step 4c
            
            BWWs = RootHelper.RefineBuds(BWWs, bShow);%%Step 4d 进一步优化芽尖。
  %             return;  

        end

        %% Step 1: 提取种子完整区域，并去除小噪声
        function BW = ExtractWholeMask(img, bShow)
% %             BW = imbinarize(img);%%  对258-1.jpg 不好，暗的芽区域太大
            imgEn = imlocalbrighten(img);%%图像局部增强
            BW = imbinarize(imgEn);
            w = 1;
            BW = Helper.DoMorphByW(BW, w);

            w = 11;%%去掉小噪点  从10增加到11
            masker = imerode(BW, strel('square', w));
            mask = BW;
            BW = imreconstruct(masker, mask);%%形态学重建
            if bShow
                figure,imshow(imoverlay(img, bwperim(BW, 8), 'red')),title('形态学重建,去噪')
                % %                 figure,imshow(BW),title('形态学重建,去噪')
            end
        end

        % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        %% Step 2a：连接 芽和种子之间的 弱连通
        function BW6 = SetConnectRegion(BW, bShow)%% for 257-1.jpg
           %% 有了Step1的imlocalbrighten，不再需要
            w = 9;
            BW2 = imdilate(BW, strel('square', w));
            BW3 = bwmorph(BW2, 'shrink', w*4);
            BW4 = BW3 & ~BW;%%获得的 额外连通分量

            CC = bwconncomp(BW);
            label = labelmatrix(CC);
            w = 4;
            BW5 = imdilate(BW4, strel('square', w));

            stats3 = regionprops('table', BW5, 'PixelIdxList');
            PixelIdxList = stats3.PixelIdxList;
            len = size(stats3, 1);
            for i = 1 : len
                pidList = PixelIdxList{i};
                labels = label(pidList);
                unLabel = unique(labels);
                if length(unLabel) <= 2
                    BW5(pidList) = false;%%%不是连通分量
                end
            end
            if false
                figure,imshow(BW5),title('BW4');
            end

            w = 4;
            BW5 = imdilate(BW5, strel('square', w));
            BW6 = BW | BW5;

            if bShow
                figure,imshow(BW6);title('BW6')
            end
% % % %             BW7 = DetailHelper.SetRGBwithBackBlack(BW6, BW);%%根的最新图 ;
% % % %             BW7 = DetailHelper.setBWColor(BW7, BW3, [255, 0, 0]);%%设置芽为蓝色
% % % %             figure,imshow(BW7);title('BW7')

        end

        %%%%%Step 2b：提取各个种子根系的单独 掩膜图像
        function BWs = GetSeedsMask(BW, bShow)
            stats = regionprops('table', BW, 'PixelIdxList', 'BoundingBox', 'Area');%%, 'SubarrayIdx', 'PixelList', 'SubarrayIdx'

            PixelIdxList = stats.PixelIdxList;
            BoundingBox = stats.BoundingBox;
            Area = stats.Area;

            mArea = median(Area);
            ins = Area > mArea/15;%%去掉过小的区域 原20

            if sum(ins) < length(ins)  %%简单优化 2022-08-01
                PixelIdxList= PixelIdxList(ins);
                BoundingBox = BoundingBox(ins, :);
                Area = Area(ins);
            end

            % % % % % % % % % % % % % % % % % % % % % % %
            boxes = ceil(BoundingBox);
            w = 3;
            boxes = boxes + [-w, -w, 2*w, 2*w];%%扩充边缘 w 像素

            lenS = length(PixelIdxList);

            BWW = cell(lenS, 1);
            for i = 1 : lenS
                pIdList = PixelIdxList{i};
                box = boxes(i, :);

                BWTmp = false(size(BW));
                BWTmp(pIdList) = true;
                BWW{i} = imcrop(BWTmp, box);
            end
            BWs.BW = BWW;
            BWs.boxes = boxes;
            BWs.lenS = lenS;

            if bShow
                figure,montage(BWW, 'Size', [1, lenS], 'BackgroundColor', 'green'),title('单独分割的种子图')
                %%'BorderSize', [1,1],
            end
        end
     
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        %% Step 3a: 截取不同颜色分量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function BWs = getColorFromBox(BWs, ccc, rrr, yyy)  %%2022-08-04增加yyy分量 来去除边缘暗分量 275-1.pg
            lenS  = BWs.lenS;
            boxes = BWs.boxes;
            BWs.ccc = cell(lenS, 1);
            BWs.rrr = cell(lenS, 1);
            for i = 1 : lenS
                box = boxes(i, :);
                BWs.ccc{i} = imcrop(ccc, box);%%这里应该会出现 其他种子的颜色分量也落入 box的情况 2022-08-01
                BWs.rrr{i} = imcrop(rrr, box);%%但后面会用 BW做掩膜！
                BWs.yyy{i} = imcrop(yyy, box);%%但后面会用 BW做掩膜！
            end
        end
         %% Step 3b New: 进行阈值分割，获得多个种子区域的不同分割成分%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         %% 原来的BWseedBud，BWroot，BWrootBud，BWseedRoot 放到 Step 3d去做  2022-04-04
         %% ps：对于259-5.jpg last种子，其种子区域 包含了较多 芽区域
         %%     对于259-4.jpg 9th种子，其种子区域 包含了较多 芽区域
         function BWs = getSeedRegion(BWs, bShow)
            
            threshC3 = BWs.threshC3;
            threshR3 = BWs.threshR3;
            threshC2 = BWs.threshC2;
            threshR2 = BWs.threshR2;
            
            lenS = BWs.lenS;
            BWs.BWseed = cell(lenS, 1);
            
            for i = 1 : lenS
                ccc = BWs.ccc{i};%%ccc=aaa
                rrr = BWs.rrr{i};%%rrr=bbb

                BW = BWs.BW{i};
                % % %  2022-08-03 老的seed区域  % % % % % %
                % % % %                 BWr = ccc >= threshC3(2);%%2022-05-24
                % % % %                 BWseed = BWr & BW;
                
                % % % 2022-08-03 改善Seed区域的搜素，从step3d改过来  % % % % % %
                BWr1 = rrr < threshR2;% 根2 是否以后采用局部的 阈值分离根效果更好！
                BWsb = ~BWr1 & BW;% 去掉Root的掩膜

                threshNoR = multithresh(ccc(BWsb), 1);% % %  还是采用局部的，改为全局计算没用，过不了4-3两行.jpg
                BWc3 = (ccc >= threshNoR) & BWsb;
                
                BWseed  = BWc3 & BW;
                % % % % % % % % % % % % % % % % % % % % %
                
                BWs.BWseed{i} = BWseed;
            end
            if bShow
                figure,montage(BWs.BWseed, 'Size', [1, lenS]),title('Seed的第一次潜在区域')               
            end
         end
        
        %%Step 3c: 对连通的两个种子区域进行分割%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %% 部分排序步骤可以不要了，在step 3f中将按照坐标重排种子区域 2022-08-03进行重构
        function BWs = RefineSeedRegion(BWs, bShow)
            lenS = BWs.lenS;

            strPath = 'F:\华农水稻\root2\TLS';
            if isfolder(strPath) == 1  %%增加 app UI程序的路版型
                try
                    addpath(strPath);
                catch
                    disp('.........')
                end
            end       
           
            mergeSeed = false;
            BWs.noSeed = zeros(lenS, 1);
            for i = 1: lenS
                BW = BWs.BW{i};
                BWseed = BWs.BWseed{i};

                w = 24;
                BWseed2 = DetailHelper.DoMorph(BWseed, w) | BWseed;%%种子区域加连通
                BWseed3 = imerode(BWseed2, strel('square', 16));

% % %                 stats = regionprops('table', BWseed3, 'PixelIdxList', 'BoundingBox', 'Area', 'Centroid');
                stats = regionprops('table', BWseed3, 'Area', 'Centroid');%%暂时不用上面的，因为后面用不到
                Area = stats.Area;
                mm = max(Area)/15;
                ins = Area >= mm;
                stats = stats(ins, :); %% 去掉过小的区域                
                
                lenSeed = height(stats);
                if lenSeed == 0
                    BWs.noSeed(i) = 1;%%说明找不到种子,后面也没用到
                    continue;
                elseif lenSeed ==1
                   continue;
                end                
               
                Centroid = round(stats.Centroid);
                %%%%%判断有没有 落在同一Seed区域的 局部子域%
                Area = stats.Area;
                LLL = bwlabel(BWseed2);
                label = zeros(lenSeed, 1);
                
                for j = 1 : lenSeed
                    cent = Centroid(j, :);%两个种子的中心位置
                    label(j) = LLL(cent(2), cent(1));
                end 
                uuu = unique(label);
                
                if length(uuu) == 1
                   continue; 
                end
                
                cents = [];
                for lll = uuu'
                    ins1 = label == lll;
                    areas = Area(ins1);
                    
                    [mm, ins2] = max(areas);
                    cent = Centroid(ins1, :);
                    cents = [cents; cent(ins2, :)];                    
                end
                Centroid = cents;
                lenSeed = size(Centroid, 1);
                Area = [];
                stats = [];%%防止后面用到
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                % % % %                 BoundingBoxS = stats.BoundingBox;%%发现没有 先注释
                % % % %                 PixelIdxListS = stats.PixelIdxList;
                if lenSeed >2
                   disp('种子区域有三个连通在一起，不能处理！')
                   return
                end
                                
                mergeSeed = true;%%说明有种子合并了    
                cent1 = Centroid(1, :);%两个种子的中心位置
                cent2 = Centroid(2, :);

                distS1 = bwdistgeodesic(BW, cent1(1), cent1(2));%% 以种子1的坐标为起点的测地距, 'cityblock'
                distS2 = bwdistgeodesic(BW, cent2(1), cent2(2));%% 以种子2的坐标为起点的测地距, 'cityblock'

                distS = distS1 + distS2;%%距离 两个种子的距离和，
                minDist = min(min(distS));
                IminDist = (distS <= (minDist +1) );%%如果增加 1，下面的w也要增加


                corners = detectHarrisFeatures(IminDist, 'MinQuality', 0.1, 'FilterSize', 15);%%要21以上，保证第4角点比第4角点要显著分别
                metric = corners.Metric;
                Big = maxk(metric, 2);
                if length(Big) == 2
                   disp(['头两个角点的取值：', num2str(Big(1)/Big(2))] )%%用作检查角点的 可区分性
                end

                if bShow
                    EEE = DetailHelper.SetRGBwithBackBlack(IminDist, BW);%%根的最新图 ;
                    EEE = DetailHelper.setBWColor(EEE, BWseed3, [0, 0, 255]);%%设置芽为蓝色
                    figure,imshow(EEE),title('距离图加4角点')
                    hold on;
                    plot(corners.selectStrongest(4));
                    strong = corners.selectStrongest(1);
                    [strong.Location(1), strong.Location(2)]
                    plot(strong.Location(1), strong.Location(2), 'or', 'MarkerSize', 15);
                    hold off
                end
              
                strongLoc = corners.selectStrongest(1).Location;%%最强角点 认为是两个芽的相交处
                w = 7;
                
% % % %                 if bShow
% % % %                     circ = drawcircle('Center',strongLoc,'Radius',w);%%从角点创造圆形区域
% % % % % % % % %                     BWmask = createMask(circ) & BW; %%得到圆形区域的小区域
% % % %                 end   
                
                BWmask = DetailHelper.CreateCirclesMask(size(BW), round(strongLoc), w) & BW;%%得到圆形区域的小区域
                
                
                BWrest = BW & ~IminDist;%%不包含连通路径的剩余部分
                CC  = bwconncomp(BWrest);%%获得连通分量
                CC2 = bwconncomp(IminDist & ~BWmask);%%不连通的 两个分支
                CC3 = bwconncomp(BWmask);%%圆形区域的连通区域
                %%三区域合并
                CC.NumObjects = CC.NumObjects + CC2.NumObjects + CC3.NumObjects;
                CC.PixelIdxList = cat(2, CC.PixelIdxList, CC2.PixelIdxList, CC3.PixelIdxList);

                L = labelmatrix(CC);%%组成颜色

                lenComp = CC.NumObjects;
                RGB = label2rgb(L);
                if bShow
                    figure,imshow(RGB),title('不同颜色的连通区域')
                end

                set1 = CC2.PixelIdxList{1};
                set2 = CC2.PixelIdxList{2};

                [BW1, mask1, masker1] = RootHelper.getMaskAndMasker(set2, set1, BW);%%形态学重构
                [BW2, mask2, masker2] = RootHelper.getMaskAndMasker(set1, set2, BW);

                if bShow
                    figure,montage({BW1, BW2}, 'Size', [1, 2]),title('两个分量的分离')
                end
                BW3  = BW1 & BW2;%%两部分分支 重构后相交的地方
                BW11 = BW1 & ~BW2;%%得到不相交的部分
                BW22 = BW2 & ~BW1;

                % % % % % % % % % % % % % % % % % % % % % %
% %                 BB = bwskel(BW3 | BWmask);
                BBskel = bwskel(BW3 | IminDist);%%两个分支 + 相交部分， 得到骨架终点、分支点
                figure,imshow(BBskel),title('骨架')
                CC = bwmorph(BBskel, 'endpoints');%%骨架上的全部终点
                [rows, cols] = find(CC);
                if bShow
                    hold on
                    plot(cols, rows, 'co', 'MarkerSize', 15, 'LineWidth', 3)
                    drawcircle('Center',strongLoc,'Radius',w);%%从芽相交处(角点)创造圆形区域
                    hold off
                end
                
                values = BWseed3(CC);
                ins = values == 0;%%找出不在种子区域上的两个终点
                conLoc = [cols(ins), rows(ins)];%%牙尖终点
                seeLoc = [cols(~ins), rows(~ins)];%%种子区域 终点

                BBbranch = bwmorph(BBskel, 'branchpoints');%%全部  分支点
                [rowBra, colBra] = find(BBbranch);
                
                if bShow
                    hold on
                    plot(colBra, rowBra, 'rs', 'MarkerSize', 10)
                    hold off
                end
                
                inddd = sub2ind(size(BBbranch), rowBra, colBra);
                distBB = distS(inddd);%%快方法
% % %                 distBB = distS(BBbranch);%%慢方法
                [~, ii] = max(distBB);
                BranLocMax = [colBra(ii), rowBra(ii)];
                [~, jj] = min(distBB);
                BranLocMin = [colBra(jj), rowBra(jj)];%%最近的点

                [BWxdata1, rowy1, colx1, ~] = RootHelper.getHalfBranch(BBskel, seeLoc(1, :), BranLocMin);
                [BWxdata2, rowy2, colx2, ~] = RootHelper.getHalfBranch(BBskel, seeLoc(2, :), BranLocMin);

                if bShow
                    hold on
                    plot(colx1, rowy1, 'yo', 'MarkerSize', 2)%%画出两个分支点
                    plot(colx2, rowy2, 'mo', 'MarkerSize', 2)
                    hold off
                end

                [~, P1] = fit_2D_data(colx1, rowy1, 'no');%%TLS 计算TLS
                [~, P2] = fit_2D_data(colx2, rowy2, 'no');
                Yhat1 = polyval(P1, conLoc(:, 1));%%计算两个顶点的误差
                Yhat2 = polyval(P2, conLoc(:, 1));
                disp(Yhat1 - conLoc(:, 2))
                disp(Yhat2 - conLoc(:, 2))

                if bShow
                    hold on
                    plot(conLoc(:, 1), Yhat1, 'r+', 'MarkerSize', 8)%%画出估计值
                    plot(conLoc(:, 1), Yhat2, 'r*', 'MarkerSize', 8)%%画出估计值
                    hold off
                end
                
                hats = [(Yhat1 - conLoc(:, 2)).^2 ; (Yhat2 - conLoc(:, 2)).^2];
                mm = [1, 0, 0, 1;  0, 1, 1, 0];%%不是2维的就要报错
                errs = mm * hats;%%计算误差

                [~, in] = sort(errs);
                conLoc = conLoc(in, :);%%调整 角点次序

                [newBW1, insG1] = RootHelper.getNewBranch(BW1, conLoc(1, :), BW11);
                [newBW2, insG2] = RootHelper.getNewBranch(BW2, conLoc(2, :), BW22);

                if bShow
                    figure,montage({insG1, insG2}, 'Size', [1, 2]),title('两个顶点分量的测地距')
                    figure,montage({newBW1, newBW2}, 'Size', [1, 2]),title('最终分支结果')
                end                
                
                BWs.BW{i}     = newBW1;%%2022-08-03 不再需要保持顺序，step3f会处理
                BWs.BW{end+1} = newBW2;
                BWs.ccc{end+1}  = BWs.ccc{i};
                BWs.rrr{end+1}  = BWs.rrr{i};
                BWs.yyy{end+1}  = BWs.yyy{i};
                BWs.boxes(end+1, :) = BWs.boxes(i, :);
            end
            
            BWs.lenS = length(BWs.BW);
            BWs.BWseed = [];%%清空BWseed，需要重新找种子
        end
        
        %%%通过分支与测地距 获得完整扩展的分支区域
        function [newBW, insG] = getNewBranch(BWbranch, conLoc, BWsep)
            distF = bwdistgeodesic(BWbranch, conLoc(1), conLoc(2));%% 利用测地距, 'cityblock'
            distE = bwdistgeodesic(BWbranch, BWsep);%% 利用测地距, 'cityblock'
            distG = distE + distF;
            minG = min(min(distG));
            w = 6;
            insG = (distG <= minG + w);
            newBW = insG | BWsep;
        end


        %% 利用测地距计算各个分支与种子距离，取一半长度用作后续回归
        function [BW2, rowy, colx, distSS] = getHalfBranch(BBskel, seeLoc, BranLocMin)

            distSS = bwdistgeodesic(BBskel, seeLoc(1), seeLoc(2));%% 利用测地距, 'cityblock'
            distRM = distSS(BranLocMin(2), BranLocMin(1));

            BW2 = (distSS>=distRM/2) & (distSS<=distRM);
            [rowy, colx] = find(BW2);
        end


        %% 获取 分支
        function [BW, mask, masker] = getMaskAndMasker(setF, setT, BW)
            mask = BW;
            mask(setF) = false;
            masker = false(size(BW));
            masker(setT) = true;
            BW = imreconstruct(masker, mask);
        end
        
        %%Step 3d：提出 芽与根的区域  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function BWs = getBudAndRoot(BWs, bShow, args)
            lenS = BWs.lenS;
            threshC3 = BWs.threshC3;
            threshR3 = BWs.threshR3;
            threshC2 = BWs.threshC2;
            threshR2 = BWs.threshR2;
% %             threshNoR = BWs.threshNoR;%%全局阈值没用

            BWs.BWseed = cell(lenS, 1);
            BWs.BWseedBud = cell(lenS, 1);
            BWs.BWroot = cell(lenS, 1);
            BWs.BWrootBud = cell(lenS, 1);
            BWs.BWseedRoot = cell(lenS, 1);
            BWs.BWbud = cell(lenS, 1); %%这个没用 后面会删掉
            
            BBB = cell(lenS, 1);
            CCC = cell(lenS, 1);
            DDD = cell(lenS, 1);
            for i = 1 : lenS
                ccc = BWs.ccc{i};
                rrr = BWs.rrr{i};
% %                 yyy = BWs.yyy{i};

                BW = BWs.BW{i};                   

% % % % 2022-05-24 采用二阈值有问题，后采用三阈值                  
% % %                 BWc = ccc < threshC;
% % %                 BWr = rrr < threshR;
% % %                 BWseed = ~BWr & BW; %还是重新来吧
% % %                 BWrootBud = BWr & BW;
% % %                 BWseedBud = BWc & BW;
% % %                 BWroot = ~BWc & BW;
% % %                 BWseedRoot = BWroot | BWseed;%%%这里不能直接用，因为seed要close/open多一点

% % % % % % 三阈值挺好的，但对 7.2-4-3第二行的图片识别出处Seed，2022-08-03
% % % % % %                 BWc3 = ccc > threshC3(2);% 种-黄色，
% % % % % %                 BWc1 = ccc < threshC3(1);% 芽-绿色，
% % % % % %                 BWc2 = ccc >= threshC3(1) & ccc <= threshC3(2);% 根
% % % % % %                 
% % % % % %                 BWr1 = rrr < threshR2;% 根        

% % % % % 新尝试 2022-08-03
                BWc2 = ccc >= threshC3(1) & ccc <= threshC3(2);% 根1
                BWr1 = rrr < threshR2;% 根2 是否以后采用局部的 阈值分离根效果更好！
                
                BWsb = ~BWr1 & BW;% 去掉Root的掩膜
                
% % %                 DetailHelper.showHistgram(ccc, BWsb, '部分直方图分割 ccc', true, args.fileSaveCCCPPP)                

                threshNoR = multithresh(ccc(BWsb), 1);% % %  还是采用局部的，改为全局计算没用，过不了4-3两行.jpg
                BWc1 = (ccc < threshNoR) & BWsb;
                BWc3 = (ccc >= threshNoR) & BWsb;
% % % % % % % % % % % % % % % % % % % % % % % % %
                
                BWseed  = BWc3 & BW;
                BWbud   = BWc1 & BW;
                BWroot1 = BWc2 & BW;
                                
                BWroot2 = BWr1 & BW;
                BWroot  = BWroot1 & BWroot2;
                
                BWrootBud  = BWroot | BWbud;
                BWseedBud  = ~BWr1 & BW;
                
                BWseedRoot = BWseed | BWroot;%%%这里不能直接用，因为seed要close/open多一点
                
                BWs.BWseed{i}     = BWseed;%还是重新来吧
                BWs.BWseedBud{i}  = BWseedBud;
                BWs.BWroot{i}     = BWroot;
                BWs.BWrootBud{i}  = BWrootBud;
                BWs.BWseedRoot{i} = BWseedRoot;
                BWs.BWbud{i} = BWbud;%%2022-08-01 留待后用
                
% % %                 BBB{i} = BWseedBud & ~BWrootBud;
% % %                 CCC{i} = BWseed & ~BWroot;
% % % % % %                 BBB{i} = BWroot1;
% % % % % %                 CCC{i} = BWroot2;
% % % % % %                 DDD{i} = BWroot1 & BWroot2;
            end
            if bShow
                figure,tiledlayout(2,2,'TileSpacing','none','Padding', 'none');
                nexttile;montage(BWs.BWseed, 'Size', [1, lenS]),title('Seed')
                nexttile;montage(BWs.BWrootBud, 'Size', [1, lenS]),title('Root Bud')
                nexttile;montage(BWs.BWseedBud, 'Size', [1, lenS]),title('Seed Bud')
                nexttile;montage(BWs.BWroot, 'Size', [1, lenS]),title('Root')               
                
% % % %                 figure,tiledlayout(3,1,'TileSpacing','none','Padding', 'none');
% % % %                 nexttile;montage(BBB, 'Size', [1, lenS]),title('BWroot1')
% % % %                 nexttile;montage(CCC, 'Size', [1, lenS]),title('BWroot2')     
% % % %                 nexttile;montage(DDD, 'Size', [1, lenS]),title('BWroot12')
           end
        end
                
        %%step 3e:清除残留(不完整)的根与芽 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2022-08-01%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function BWs = refineFalseSeed(BWs, args)
            sumPixels = RootHelper.calcPixelsEachRegion(BWs, args);            
             
            ssum = sum(sumPixels, 2);%%行最大值
            ratios = sumPixels ./ ssum;
            rrr = max(ratios, [], 2);
            ins = rrr >= 0.9;
            
            %%需要重组所有cells
            if sum(ins) > 0
                ind = ~ins;
                BWs.BWseed     = BWs.BWseed(ind);
                BWs.BWseedBud  = BWs.BWseedBud(ind);
                BWs.BWroot     = BWs.BWroot(ind);
                BWs.BWrootBud  = BWs.BWrootBud(ind);
                BWs.BWseedRoot = BWs.BWseedRoot(ind);
                BWs.BWbud      = [];%%没用清空
                
                BWs.BW = BWs.BW(ind);                
                BWs.boxes = BWs.boxes(ind, :);
                BWs.ccc = BWs.ccc(ind);
                BWs.rrr = BWs.rrr(ind);
                BWs.yyy = BWs.yyy(ind);
                BWs.lenS = length(BWs.BW);
            end
        end
        
        function sumPixels = calcPixelsEachRegion(BWs, args)
            lenS = BWs.lenS;
            sumPixels = zeros(lenS, 3);
            for i = 1 : lenS                
                BWseed = BWs.BWseed{i};
                BWroot = BWs.BWroot{i};
                BWbud = BWs.BWbud{i};
                
                sumPixels(i, 1) = sum(sum(BWseed));
                sumPixels(i, 2) = sum(sum(BWroot));
                sumPixels(i, 3) = sum(sum(BWbud));                
            end         
            if args.bSaveRatio == true                
                save(args.fileSaveRatio, 'sumPixels', '-ascii', '-tabs');                
            end
        end
        
        %%step 3f: 对各个种子区域进行位置重排序 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2022-08-02%%%%%%%%%%%%%%%%%
        function BWs = reNumberSeeds(BWs, bShow, args)
            lenS = BWs.lenS;

            boxes = BWs.boxes;
            majLength = boxes(:, 4);
            centroid  = [boxes(:, 1) + 0.5 * boxes(:, 3), boxes(:, 2) + 0.5 * boxes(:, 4)];
            
            stats.Centroid = centroid;
            stats.MajorAxisLength = majLength;
            
            [ind, ~, ~, bSort] = DetailHelper.findIndexForSeed(stats);%%重置顺序
            
            if bSort == true
                BWs.BWseed     = BWs.BWseed(ind);
                BWs.BWseedBud  = BWs.BWseedBud(ind);
                BWs.BWroot     = BWs.BWroot(ind);
                BWs.BWrootBud  = BWs.BWrootBud(ind);
                BWs.BWseedRoot = BWs.BWseedRoot(ind);  

                BWs.BW = BWs.BW(ind);                
                BWs.boxes = BWs.boxes(ind, :);
                BWs.ccc = BWs.ccc(ind);
                BWs.rrr = BWs.rrr(ind);
                BWs.yyy = BWs.yyy(ind);
            end
            
            if bShow && args.bSaveNumber
                centroid = centroid(ind, :);
                figure,imshow(args.img8);
                oneone = ones(size(centroid));
                
                ptext = centroid + oneone .* [1.3, -2] *16; %%加上偏移量
                i = 1 : lenS;
                ii = num2str(i');
                text(ptext(:, 1), ptext(:, 2), ii, 'FontSize',25, 'Color', [0.4, 0.9, 0.6])%%画种子序号
                
                if args.bSaveNumber
                    DetailHelper.SaveFigure2Img(args.fileSaveNumber);
                end
            end
        end
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % Step 4a 获得种子的 扩大区域，以及 芽的初步区域 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function BWs = getBudAndSeed(BWs, bShow)
            lenS = BWs.lenS;
            BWBudss = cell(lenS, 1);

            BWseedNew = cell(lenS, 1);
            BBB = cell(lenS, 1);
            CCC = cell(lenS, 1);
            DDD = cell(lenS, 1);
            EEE = cell(lenS, 1);
            for i = 1 : lenS
                BWseed = BWs.BWseed{i};
                BWseedBud = BWs.BWseedBud{i};
                BWroot = BWs.BWroot{i};
                BW = BWs.BW{i};

                % % % % % 先优化root 2022-08-04 % % % % % % % % % % % % % % % % % %
                yyy = BWs.yyy{i};
                ddd = yyy(BW);
                mmm = mean(ddd);
                sss = std(ddd);
                thresh = mmm - sss;
                BWr = yyy > thresh;
                BWroot = BWroot & BWr;
                BWs.BWroot{i} = BWroot;                
                % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
                
                %%% 第三种做法：
                w  = 2;
                BWseed2 = imdilate(BWseed, strel('square', w)) & BW;
                BWroot2 = imdilate(BWroot, strel('square', w)) & BW;
                w = 16;
                BWseed2 = DetailHelper.DoMorph(BWseed2, w) | BWseed;%%种子区域加连通
% %                 BWseed2 = imclose(BWseed2, strel('square', w));%%种子区域加连通
%                 DetailHelper.AreaFilter(BWseed2, BW)

                BWseedRoot = BWseed2 | BWroot2;%%种子与根区域
                DDD{i} = BWseedBud & ~BWseedRoot;%%留下芽

                BWseedNew{i} = BWseed2;%%保存最新的Seed
                
% %                 figure,montage({BW & ~BWseedBud,BWroot , DDD{i} }, 'Size', [1, 3]),title('对比BW与Refined芽的差异')
                BWBudss{i} = DDD{i};
                
                BBB{i} = BWseedRoot;
               
            end
            if bShow
%                 figure,montage(BWBudss, 'Size', [1, lenS]),title('最新的芽图')
                figure,montage(BBB, 'Size', [1, lenS]),title('种子所在位图'),DetailHelper.SetWindowTitle();
                figure,montage(BWseedNew, 'Size', [1, lenS]),title('Seed图2'),DetailHelper.SetWindowTitle();
                figure,montage(DDD, 'Size', [1, lenS]),title('芽图3'),DetailHelper.SetWindowTitle();
            end

            BWs.BWseedNew = BWseedNew;
            BWs.BWbudNew = BWBudss;
        end

        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %  
        %%% Step 4b：完善芽的区域  % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
        function BWs = getRefinedBuds(BWs, bShow)
            lenS = BWs.lenS;

            BWBudss = cell(lenS, 1);
            BBB = cell(lenS, 1);
            CCC = cell(lenS, 1);
            DDD = cell(lenS, 1);
            EEE = cell(lenS, 1);
            FFF = cell(lenS, 1);
            for i = 1 : lenS
                BW = BWs.BW{i};
                BWbud = BWs.BWbudNew{i};
                BWseed = BWs.BWseedNew{i};

                % % % % % % % % % % % % % % % % Step 1，通过测地距找出连通区域
                BWseed = bwareafilt(BWseed, 1, 8);%%取最大的连通区域  %%%2022-05-24 去最大的两个连通区域求凸包
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
                BWs.BWseedNew{i} = BWseed;%%重新设置最新的seed区域！！

                BWbud2 = bwareafilt(BWbud, [15, Inf]);%%去掉较小的连通区域
                BWbud2 = imclose(BWbud2, strel('square', 8));%%链接一些分离的区域      

                [row, col] = find(BWseed);
                distS = bwdistgeodesic(BW, col, row, 'cityblock');%% 以种子区域为起点的测地距, 增大种子的 范围

                stats = regionprops('table', BWbud2, 'Area', 'PixelList', 'PixelIdxList');%%%, 'PixelList'
                PixelList = stats.PixelList;
                PixelIdxList = stats.PixelIdxList;
                Area = stats.Area;
                distError = 8;%%%容忍的宽度

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%最新版本：for 257-1.jpg，因为芽到根有些地方会 不连通
                BWbud3 = BWbud2;
                
                for jj = 1 : length(Area)
                    pList   = PixelList{jj};
                    pIdList = PixelIdxList{jj};

                    distB = bwdistgeodesic(BW, pList(:, 1), pList(:, 2), 'cityblock');%% 利用测地距, 增大芽的 范围
                    minBudToSeed = min(distS(pIdList));%%计算最短的测地距

                    dist = distS + distB;
 
                    BWins1 = dist <= (minBudToSeed + 1); 
                    BWins2 = BWins1 & (dist >=8);
                    BWins3 = imdilate(BWins2, strel('square', 8));
                    BWins  = BWins1 | BWins3;
                    
                    BWbud3 = BWbud3 | BWins;%%得到芽区域
                end
                BWbud3 = BWbud3 & ~BWseed;
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 连通芽种的两主干
                stats = regionprops('table', BWbud3, 'Area', 'PixelList', 'PixelIdxList');%%%, 'PixelList'
                if height(stats)>1 && 1 == 1
                    Area = stats.Area;
                    
                    [~, ins] = sort(Area,'descend');
                    stats = stats(ins, :);
                    
                    PixelList = stats.PixelList;
                    PixelIdxList = stats.PixelIdxList;
                    Area = stats.Area;
                    
                    maxA = Area(1);
                    
                    pListMax = PixelList{1};
                    distM = bwdistgeodesic(BW, pListMax(:, 1), pListMax(:, 2), 'cityblock');%% 利用测地距,
                    for jj = 2 : length(Area)
                        pList   = PixelList{jj};
                        pIdList = PixelIdxList{jj};
                        
                        distB = bwdistgeodesic(BW, pList(:, 1), pList(:, 2), 'cityblock');%% 利用测地距, 增大芽的 范围
                        
                        minBudToSeed = min(distM(pIdList));%%计算最短的测地距
                        
                        dist = distM + distB;
                        
                        BWins1 = dist <= (minBudToSeed + 1);
                        BWbud3 = BWbud3 | BWins1;%%得到芽区域
                    end
                end                
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                  
                BWbud4 = BWbud3;
                CCC{i} = imoverlay(BWbud4, bwperim(BWseed, 8), 'yellow');

                % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
                % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
                % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
               %% Step 3:扩展空白区域，不OK：258-1.jpg，下面大改
% % % % % %                 BWseedRoot = BWs.BWseedRoot{i} | BWseed;  %%去掉散根，好和芽合并
% % % % % % % % % % % %                 BWseedRoot = imclose(BWseedRoot, strel('square', 3));
% % % % % %                 masker = imerode(BWseedRoot, strel('square', 8));
% % % % % %                 BWseedRoot = imreconstruct(masker, BWseedRoot);
% % % % % %                 mask = BW & ~BWseedRoot;
% % % % % %
% % % % % %                 BWbud6 = imreconstruct(BWbud4, mask) | BWbud4;                EEE{i} = BWbud6;
% % % % % %
% % % % % % % % % % % %                 figure,montage({BW & ~BWbud4,BW & ~BWbud6, BWbud6}, 'Size', [1, 3]),title('对比BW与Refined芽的差异2')
% % % % % %                 BWBudss{i} = BWbud6;
% % % % % %
% % % % % %                 FFF{i} = DetailHelper.SetRGBwithBackBlack(BWbud6, BW);
                % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
                % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
                % % % % % % % Step 3：OK, for 258-1.jpg但一些底部的宽度不够，很多芽的尖尖处还没充满， 需要得到“得到长芽部分”
                [row, col] = find(BWbud4);
                distB = bwdistgeodesic(BW, col, row, 'cityblock');%%计算与芽的测地距
                
               %% {2022-05-22： OK但不好写论文，需要进一步改
                distC = distS - distB;
                BWcontour = imdilate(BWbud4, strel('square', 2)) & ~BWbud4;%%扩大轮廓
                distC2 = rmmissing(distC(BWcontour));%%去除NaN的值
                meanC2 = median(distC2);%%取中值
                ins = distC > meanC2;%maxSeed * 0.95;                
                %%}

                BWbud6 = BWbud4 | ins;
                EEE{i} = BWbud6;
                BWBudss{i} = BWbud6;    
                    
            end
            if bShow
% %                 figure,montage(BBB, 'Size', [1, lenS]),title('最新的芽图')
                figure,montage(CCC, 'Size', [1, lenS]),title('芽图-BWbud4'),DetailHelper.SetWindowTitle();
% %                 figure,montage(DDD, 'Size', [1, lenS]),title('最新的芽图44'),DetailHelper.SetWindowTitle();
                figure,montage(EEE, 'Size', [1, lenS]),title('芽图455-BWbud6'),DetailHelper.SetWindowTitle();
% %                 figure,montage(FFF, 'Size', [1, lenS]),title('最新的芽图-彩色'),DetailHelper.SetWindowTitle();
            end
            BWs.BWbudNew = BWBudss;

        end
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
        %%% Step 4c: 提取主根 以及完整的 主根区域 % % % % % %  % % % % % % % % % % % % % % % % % % 
        function BWs = getMajorRoots(BWs, bShow)
            lenS = BWs.lenS;

            BWRootss = cell(lenS, 1);%% 保存主根
            originRoot = zeros(lenS, 2);%%保存芽种的交界
            BWbud2root = cell(lenS, 1);%%保存芽 以及与种的交界处
            
            BBB = cell(lenS, 1);
            CCC = cell(lenS, 1);
            DDD = cell(lenS, 1);
            EEE = cell(lenS, 1);
            for i = 1 : lenS
                BW = BWs.BW{i};
                BWbud  = BWs.BWbudNew{i};
                BWseed = BWs.BWseedNew{i};
                %                 BWseedBud = BWs.BWseedBud{i};
                BWroot = BWs.BWroot{i};
                %                 BWrootBud = BWs.BWrootBud{i};
                %                 BWseedRoot = BWs.BWseedRoot{i};
                ccc = BWs.ccc{i};

                % % % % % % % % % % % % % % % % % % % % % % % %
                %%%% Step 1：去除弱根
                mask = BWroot & ~BWbud;
                mask = imclose(mask, strel('square', 4)) | mask;
                masker = imerode(mask, strel('square', 7));

                BWroot2 = imreconstruct(masker, mask);
% %                 figure,montage({BW & ~BWroot2, BWroot, BWroot2}, 'Size', [1, 3]),title('对比BW与根的差异')

                BBB{i} = DetailHelper.SetRGBwithBackBlack(BWroot2, BW);%%根的最新图


                % % % % % % % % % % % % % % % % % % % % % % % %
                %%%%Step 2：构建主根，如果主根不靠近 芽种边界，则采用graydist来求出路径，并在种子上找出弱根！
                w = 4;
                BWbudEx  = imdilate(BWbud, strel('square', w)) & BW;
                BWseedEx = imdilate(BWseed, strel('square', w)) & BW;
                BWborderSeedBud =  BWbudEx & BWseedEx;%%构建 芽种边界！
                
%                 BWbud2root{i} = BWbud | BWborderSeedBud; %%保存芽，以及交界处

% % % % % %                 CCC{i} = DetailHelper.setBWColor(BBB{i}, BWbud, [0, 0, 255]);%%设置芽为蓝色
% % % % % %                 CCC{i} = DetailHelper.setBWColor(CCC{i}, BWborderSeedBud, [255, 0, 0]);%%设置芽种边界为红色

                % % % % % % % % % % % % % % % % % % % % % % % % %
                %% 因为 55-2.jpg 有两条主根差别不大，需要额外处理，找出主根
                stats = regionprops('table', BWroot2, 'Centroid', 'PixelIdxList', 'Area');
                if height(stats) >1 %% 去掉假根，距离芽最近的不是主根
                    area = stats.Area;
                    [B, in] = maxk(area, 2);
                    
                    if B(2) > B(1) *2/3
                        stats = stats(in, :);
                        area = stats.Area;
                        cent = stats.Centroid;
                        pList = stats.PixelIdxList;
                        
                        stats = regionprops('table', BWbud, 'centroid');
                        if height(stats) > 1
                            disp('error in regionprops of BWbud')
                        end
                        centBud = stats.Centroid;
                        
                        dist1 = norm(cent(1, :) - centBud);
                        dist2 = norm(cent(2, :) - centBud);
                        
                        if dist2 > dist1 * 2
                            BWroot2(pList{1}) = false; %% 去掉最大的根
                            disp('erase the large false root')
                        end                        
                    end
                end
                
                % % % % % % % % % % % % % % % % %
                BWrootMajor = bwareafilt(BWroot2, 1);%%主根区域

                meanSeed = mean2(ccc(BWroot2));%%获取根部的 统计特征
                ccc2 = (ccc - meanSeed) .^2;%%平方增大差异
                ccc2(~BW) = Inf;

% % % % % % % % % % % % % % % % % % % % % %
                BWshink = bwmorph(BWborderSeedBud, 'shrink');%%形态学，将收缩到中点
                [rowBud, colBud] = find(BWshink);  %%OK
                meanRow = round(mean(rowBud));
                meanCol = round(mean(colBud));

                T2 = graydist(ccc2, meanCol, meanRow);%%灰度距离，从中点出发
% % % % %                 T2 = graydist(ccc2, BWborderSeedBud);%%灰度距离OK， 被上面替代
                T1 = graydist(ccc2, BWrootMajor);%%灰度距离，从根部出发

                TT = T1 + T2;
                minGray = min(min(TT));%%最短灰度距离
                BWrootExtent = TT <= minGray * 1.03;%%得到根的扩充，可扩充到芽种边界，原1.05

                BWrootMajor = imreconstruct(BWrootMajor, BWroot2 | BWrootExtent);

                if bShow
                    EEE{i} = DetailHelper.SetRGBwithBackBlack(BWrootMajor, BW);%%根的最新图 ;
                    EEE{i} = DetailHelper.setBWColor(EEE{i}, BWbud, [0, 0, 255]);%%设置芽为蓝色
                    EEE{i} = DetailHelper.setBWColor(EEE{i}, BWborderSeedBud, [255, 0, 0]);%%设置芽种边界为红色
                    EEE{i} = DetailHelper.setBWColor(EEE{i}, BWrootExtent & ~BWroot2, [255, 10, 250]);%%设置扩充根部为紫色
                end

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
              
                BWRootss{i} = BWrootMajor; %% 保存主根
                originRoot(i, :) = [meanCol, meanRow];%%保存芽种分出口坐标
            end


            if bShow
                figure,montage(BBB, 'Size', [1, lenS]),title('根图4c-1'),DetailHelper.SetWindowTitle();
% % % % %                 figure,montage(CCC, 'Size', [1, lenS]),title('最新的根图2')
                fig = figure;montage(EEE, 'Size', [1, lenS]),title('根图4c-2'),DetailHelper.SetWindowTitle();
% %                 figure,montage(DDD, 'Size', [1, lenS]),title('最新的根图4')
                BWs.fig = fig;
            end
            BWs.BWrootMajor = BWRootss;
            BWs.originRoot = originRoot;
            BWs.BWbud2root = BWbud2root;
        end
        
        %%% Step 4d: 进一步强化芽区域，主要强化芽远端，如果不用主干计算长度，那其实用途不大 %%%%%%%%%%%%%%%%%%%%%%%
         function BWs = RefineBuds(BWs, bShow)
            lenS = BWs.lenS;
           
            BBB = cell(lenS, 1);
            CCC = cell(lenS, 1);
            DDD = cell(lenS, 1);
            EEE = cell(lenS, 1);
            
            BWbudNewEx = cell(lenS, 1);
            for i = 1 : lenS
                BW = BWs.BW{i};
                BWbud  = BWs.BWbudNew{i};
                BWseed = BWs.BWseedNew{i};
                BWrootMajor = BWs.BWrootMajor{i};

                BWseedEx = imdilate(BWseed, strel('square', 16));
                
                BWseedEx = bwareafilt(BWseedEx, 1);%%主种区域
                stats = regionprops('table', BWseedEx, 'MajorAxisLength', 'MinorAxisLength');
                MajorAxisLength = stats.MajorAxisLength;
                MinorAxisLength = stats.MinorAxisLength;
                
                %%radii = round((MajorAxisLength + MinorAxisLength) / 2) - 16;
                radii = round(MinorAxisLength - 16);
                
                if radii > 0
                    BWseedEx = imdilate(BWseedEx, strel('square', radii));
                end
                
                % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
                
                BWbudEx  = imdilate(BWbud, strel('square', 16)) & BW;
                BWbudEx2 = BWbudEx & ~BWseedEx & ~BWrootMajor & BW;
                BWbudEx2 = BWbudEx2 | BWbud;
                
                BWbudNewEx{i} = BWbudEx2 ;
                % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
                
                BBB{i} = BWseedEx & BW;
                CCC{i} = BWbudEx;
                DDD{i} = BWbudEx2;                   
            end
            
            if bShow
                figure,montage(BBB, 'Size', [1, lenS]),title('种子imdilate-1'),DetailHelper.SetWindowTitle();
                figure,montage(CCC, 'Size', [1, lenS]),title('芽的imdilate-1'),DetailHelper.SetWindowTitle();
                figure,montage(DDD, 'Size', [1, lenS]),title('最新的芽扩充-1'),DetailHelper.SetWindowTitle();
%                 BWs.fig = fig;
            end            
            BWs.BWbudNew = BWbudNewEx;
         end
        
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % %  最后一步 : 计算芽长与根长
        function [data, figBud, figRoot, figAll] = calcRootBudLength(BWs, args)
% %             args.mmPerPixel
            [ret, figBud, BWs] = RootHelper.getBudLength(BWs, args.bShow);
            budLen = ret;
            data.budPixels = budLen;%%Step 1: 计算芽的长度
            data.budLength = budLen * args.mmPerPixel /10;  %% 转成cm单位
% %             ret(:, 2) * args.mmPerPixel

            [ret, figRoot, BWs] = RootHelper.getRootLength(BWs, args.bShow);            
            rootLen = ret;
            data.rootPixels = rootLen;%%Step 2: 计算根的长度
            data.rootLength = rootLen * args.mmPerPixel /10;

% %             ret * args.mmPerPixel
            figAll = RootHelper.DrawSeedOnGraph(args.img8, BWs);
            
            if args.saveFinal == true
                DetailHelper.SaveFigure2Img(args.saveFinalResult);                
            end
        end

        %%%Step 1: 计算芽的长度
        function [data, fig, BWs] = getBudLength(BWs, bShow)
            lenS = BWs.lenS;

            BBB = cell(lenS, 1);

            tipPos = zeros(size(BWs.originRoot));            
% %             data = zeros(lenS, 3);
            data = zeros(lenS,1);
            for i = 1 : lenS
                BW = BWs.BW{i};
                BWbud = BWs.BWbudNew{i};
% %                 BWseed = BWs.BWseedNew{i};
                originRoot = BWs.originRoot(i, :);
%                 BWbud2root = BWs.BWbud2root{i};
                
                D2 = bwdistgeodesic(BW, originRoot(1), originRoot(2), 'quasi-euclidean');%% 接口点测地距   
                
                [dist2, ins] = max(D2(BWbud)); %%老方法，OK
                [rows, cols] = find(BWbud);
                row = rows(ins);
                col = cols(ins);%%得到牙尖坐标
                
                data(i) = dist2;%%老的测地距
                
                if 1 == 0   %%OK，准确的测地距，但误差反而增大，暂时不用
                    D3 = bwdistgeodesic(BW, col, row, 'quasi-euclidean');%%牙尖的测地距
                    [dist4, BWskel] = DetailHelper.GetDistForAPath(BWbud, D2, D3);%%最新计算骨干长度
                    data(i) = dist4;%%新的骨干方法
                else
                    BWskel = bwskel(BWbud);                
                end

                if bShow
                    BBB{i} = DetailHelper.SetRGBwithBackBlack(BWbud, BW);%%根的最新图 ;

                    w = 11;
                    BWmask = DetailHelper.CreateCirclesMask(size(BW), round(originRoot), w) & BW;%%得到圆形区域的小区域
                    BBB{i} = DetailHelper.setBWColor(BBB{i}, BWmask, [255, 0, 0]);%%设置芽接口为黄色

                    w = 11;
                    BWmask = DetailHelper.CreateCirclesMask(size(BW), round([col, row]), w);%%得到圆形区域的小区域-
                    BBB{i} = DetailHelper.setBWColor(BBB{i}, BWmask, [255, 0, 0]);%%设置芽接口为黄色
          
                    BBB{i} = DetailHelper.setBWColor(BBB{i}, BWskel, [255, 0, 0]);%%设置芽为蓝色
                end
                
                tipPos(i, :) = [col, row];%%牙尖位置
            end
            if bShow
                fig = figure;montage(BBB, 'Size', [1, lenS]),title('骨干芽图'),DetailHelper.SetWindowTitle();
                % % % % %                 figure,montage(CCC, 'Size', [1, lenS]),title('最新的根图2'),DetailHelper.SetWindowTitle();
                % % % % %figure,montage(EEE, 'Size', [1, lenS]),title('最新的根图3'),DetailHelper.SetWindowTitle();
                % %                 figure,montage(DDD, 'Size', [1, lenS]),title('最新的根图4'),DetailHelper.SetWindowTitle();
            else
                fig = 0;
            end
% %             BWs.BWtrucks = BWtrucks;
            BWs.tipPos = tipPos;       %%%%  
        end
        
        %% Step2：计算根的长度
        function [data, fig, BWs] = getRootLength(BWs, bShow)
            lenS = BWs.lenS;

            BBB = cell(lenS, 1);
            CCC = cell(lenS, 1);
            DDD = cell(lenS, 1);
            EEE = cell(lenS, 1);
            data = zeros(lenS, 1);
            rootEnd = zeros(size(BWs.originRoot));
            for i = 1 : lenS
                BW = BWs.BW{i};
                BWrootMajor = BWs.BWrootMajor{i};
                BWseed = BWs.BWseedNew{i};
                originRoot = BWs.originRoot(i, :);


                D2 = bwdistgeodesic(BWrootMajor, originRoot(1), originRoot(2), 'quasi-euclidean');%%

                % % % % % % % % % % % % % % % % % % % OK，老方法：
% % %                 [dist2, ins] = max(D2(BWrootMajor));
% % %                 data(i) = dist2;      
% % %                 
% % %                 [rows, cols] = find(BWrootMajor);
% % %                 ptr = [cols(ins), rows(ins)];%%根末端
                % % % % % % % % % % % % % % % % % % 新的方法                
                [dist2, ins] = max(D2, [], 'all', 'linear');        
                [row, col] = ind2sub(size(D2), ins);%%根末端   
                ptr = [col, row];
                
                % % % % % % % % % % % % % 用最新的主干来逼近medial axis 
                if 1 == 0  %%为了提速，暂时不用这个
                    D3 = bwdistgeodesic(BWrootMajor, ptr(1), ptr(2), 'quasi-euclidean');%%
                    [dist4, BWskel] = DetailHelper.GetDistForAPath(BWrootMajor, D2, D3);%%最新计算骨干长度                
                    data(i) = dist4;%%新骨干方法
                else
                    BWskel = bwskel(BWrootMajor);
                end
                
                % % % % % % % % % % % % % % % %                
                data(i) = dist2;%%老的测地距
                
                rootEnd(i, :) = ptr;
                
                if bShow == true
                    % BWskel = bwskel(BWrootMajor);
                    BBB{i} = DetailHelper.SetRGBwithBackBlack(BWrootMajor, BW);%%根的最新图 ;
                    BBB{i} = DetailHelper.setBWColor(BBB{i}, BWskel, [255, 0, 0]);%%设置芽为蓝色
                    
                    w = 11;
                    BWmask = DetailHelper.CreateCirclesMask(size(BW), round(originRoot), w) & BW;%%得到圆形区域的小区域
                    BBB{i} = DetailHelper.setBWColor(BBB{i}, BWmask, [255, 0, 0]);%%设置芽接口为黄色
                end
            end
            if bShow                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                
                fig = figure;montage(BBB, 'Size', [1, lenS]),title('骨干根图'),DetailHelper.SetWindowTitle();
                % % % % %                 figure,montage(CCC, 'Size', [1, lenS]),title('最新的根图2'),DetailHelper.SetWindowTitle();
                % % % % %figure,montage(EEE, 'Size', [1, lenS]),title('最新的根图3'),DetailHelper.SetWindowTitle();
                % %                 figure,montage(DDD, 'Size', [1, lenS]),title('最新的根图4'),DetailHelper.SetWindowTitle();
            else
                fig = 0;
            end
            BWs.rootEnd = rootEnd;       %%%%  
        end
        %% Step3：画出图中的信息
        function figAll = DrawSeedOnGraph(img, BWs)            
            lenS = BWs.lenS;
            
            figAll = figure;%,imshow(img),tile('');
            
            for i = 1 : lenS
                BWrootMajor = BWs.BWrootMajor{i};
                BWseed = BWs.BWseedNew{i};
                BWbud = BWs.BWbudNew{i};                    
       
                box = BWs.boxes(i, :);
                      
                [rows, cols] = RootHelper.getCordsBW(BWseed, box, 3);
                img = DetailHelper.setImgWithColor(img, rows, cols, [255, 64, 0]);

                [rows, cols] = RootHelper.getCordsBW(BWrootMajor, box, 2);
                img = DetailHelper.setImgWithColor(img, rows, cols, [0, 255, 0]);
                
                [rows, cols] = RootHelper.getCordsBW(BWbud, box, 2);
                img = DetailHelper.setImgWithColor(img, rows, cols, [0, 64, 255]);         
            end
            
            imshow(img)%%画出构图
            
            hold on
            originRoot = BWs.originRoot;
            tipPos = BWs.tipPos;
            rootEnd = BWs.rootEnd;                
            boxes = BWs.boxes;
            oneone = ones(size(originRoot));
            
            originRoot = originRoot + boxes(:, 1:2) - oneone;
            tipPos = tipPos + boxes(:, 1:2) - oneone;
            rootEnd = rootEnd + boxes(:, 1:2) - oneone;
            
            ww = 7;  
            plot(originRoot(:, 1), originRoot(:, 2), 'ro', 'Markersize', ww, 'MarkerFaceColor', 'c');%%画节点
            plot(tipPos(:, 1), tipPos(:, 2), 'mo', 'Markersize', ww, 'MarkerFaceColor', 'm');
            plot(rootEnd(:, 1), rootEnd(:, 2), 'yo', 'Markersize', ww, 'MarkerFaceColor', 'y');
            
            ptext = originRoot + oneone .* [1.2, -2] *16; %%加上偏移量
            i = 1 : lenS;
            ii = num2str(i');
            text(ptext(:, 1), ptext(:, 2), ii, 'FontSize',24, 'Color', [0.4, 0.9, 0.6])%%画种子序号

            hold off
        end
        %% 根据box偏移得到BW中的有用坐标
        function [rows, cols] = getCordsBW(BW, box, ww)
            leftCol = box(1);
            leftRow = box(2);
            BW2 = bwperim(BW, 8);
            BW2 = imdilate(BW2, strel('square', ww));%%加粗边界显示
            [rows, cols] = find(BW2);
            rows = rows + leftRow - 1;
            cols = cols + leftCol - 1;  
        end
    end
end